
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
import cartopy as cp
import cartopy.feature as cfeature
import cartopy.crs as ccrs
from pylab import *


may = xr.open_dataset('../monsoon/datasets/skt.may.nc')


# may = may.skt
clim_may = may.mean(dim='time')


from cartopy.util import add_cyclic_point

clim_may_data = clim_may['skt']
lon = clim_may.coords['lon']

print("Original shape -", clim_may_data.shape)
lon_idx = clim_may_data.dims.index('lon')
wrap_clim_may, wrap_lon = add_cyclic_point(clim_may_data.values, coord=lon, axis=lon_idx)
print("New shape -", wrap_clim_may.shape)

# Plotting the data 

# lon = clim_may.lon
lat = clim_may.lat

fig = plt.figure(figsize=(16,16))
ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=0.0, globe=None))

mp = ax.imshow(wrap_clim_may-273.5,extent=(wrap_lon.min(),wrap_lon.max(),lat.min(),lat.max()),cmap='jet',origin='upper')
plt.title('Long Term Mean Temperature (May)')
# plt.legend(['Temp'])

states_provinces = cfeature.NaturalEarthFeature(
        category='cultural',
        name='admin_1_states_provinces_lines',
        scale='10m',
        facecolor='none')
# ax.add_feature(cfeature.BORDERS,edgecolor='blue')
# ax.add_feature(states_provinces, edgecolor='blue')

ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.OCEAN)


cbar = fig.colorbar(mp, shrink=0.3,label='Temperature (degC)')
cbar.minorticks_on()

#adding the long lat grids and enabling the tick labels
gl = ax.gridlines(draw_labels=True,alpha=0.1)
gl.top_labels = False
gl.right_labels = False

Original shape - (94, 192)
New shape - (94, 193)


from cartopy.util import add_cyclic_point


months = ['may','june','july','august','september','jjas']
# lon = clim_may.lon
lat = clim_may.lat
# plotting all the clim_data graphs at once

for i in months:
    clim_data = xr.open_dataset('../monsoon/datasets/skt.'+i+'.nc')
    clim_data = clim_data.mean(dim='time')

    # adding the cyclic points to remove the missing values in the center of the plot
    data = clim_data['skt']
    lon = clim_data.coords['lon']

    print("Original shape -", data.shape)
    lon_idx = data.dims.index('lon')
    clim_wrap_data, wrap_lon = add_cyclic_point(data.values, coord=lon, axis=lon_idx)
    print("New shape -", clim_wrap_data.shape)

    # clim_data = clim_data.skt
    # clim_data = clim_data.mean(dim='time')
    fig = plt.figure(figsize=(16,16))
    ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=0.0, globe=None))
    mp = ax.imshow(clim_wrap_data-273.5,extent=(wrap_lon.min(),wrap_lon.max(),lat.min(),lat.max()),cmap='jet',origin='upper')
    plt.title('Mean Surface Temperature ('+i.capitalize()+')')
    ax.add_feature(cfeature.LAND)
    ax.add_feature(cfeature.COASTLINE)
    ax.add_feature(cfeature.OCEAN)
    cbar = fig.colorbar(mp, shrink=0.3,label='Temperature (degC)')
    cbar.minorticks_on()
    gl = ax.gridlines(draw_labels=True,alpha=0.1)
    gl.top_labels = False
    gl.right_labels = False
    plt.savefig('../monsoon/figures/Surface Temperature ('+i+').png')

Original shape - (94, 192)
New shape - (94, 193)
Original shape - (94, 192)
New shape - (94, 193)
Original shape - (94, 192)
New shape - (94, 193)
Original shape - (94, 192)
New shape - (94, 193)
Original shape - (94, 192)
New shape - (94, 193)
Original shape - (94, 192)
New shape - (94, 193)


# function to add cyclic points 

def add_cyclic(data, lon):
    """Add a cyclic point to a data array.
    Parameters
    ----------
    data : array-like
        An n-dimensional array of data to add a cyclic point to.
    lon : array-like
        An array which specifies the coordinate values for
        `data`'s longitude coordinate.
    Returns
    -------
    The data array with a cyclic point added.
    """
    cyclic_data, cyclic_lon = add_cyclic_point(data, coord=lon, axis=-1)
    return cyclic_data, cyclic_lon


# Doing the same operation as above for MSLP dataset 

months = ['may','june','july','august','september','jjas']
# lon = clim_may.lon
lat = clim_may.lat
# plotting all the clim_data graphs at once

for i in months:
    clim_data = xr.open_dataset('../monsoon/datasets/mslp.'+i+'.nc')

    # adding the cyclic points to remove the missing values in the center of the plot using the above defined function
    clim_data = clim_data.mean(dim='time')

    cyclic_data, cyclic_lon = add_cyclic_point(clim_data['mslp'].values, coord=clim_data.coords['lon'], axis=-1)



    # clim_data = clim_data.mslp
    # clim_data = clim_data.mean(dim='time')
    fig = plt.figure(figsize=(16,16))
    ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=0.0, globe=None))
    mp = ax.imshow(cyclic_data/100,extent=(cyclic_lon.min(),cyclic_lon.max(),lat.min(),lat.max()),cmap='jet',origin='upper')
    plt.title('Long Term Mean MSLP ('+i.capitalize()+')')
    ax.add_feature(cfeature.LAND)
    ax.add_feature(cfeature.COASTLINE)
    ax.add_feature(cfeature.OCEAN)
    cbar = fig.colorbar(mp, shrink=0.3,label='MSLP (hPa)')
    cbar.minorticks_on()
    gl = ax.gridlines(draw_labels=True,alpha=0.1)
    gl.top_labels = False
    gl.right_labels = False
    plt.savefig('../monsoon/figures/MSLP ('+i+').png')


# Plotting the time series of skin temp and MSLP for a particular location

clim_skt = xr.open_dataset('../monsoon/datasets/skt.ymonmean.nc')
clim_mslp = xr.open_dataset('../monsoon/datasets/mslp.ymonmean.nc')


clim_mslp = clim_mslp.mslp
clim_skt = clim_skt.skt

# Plotting the time series of skin temp and MSLP for a particular location

clim_skt_land = clim_skt.sel(lat=20,lon=80,method='nearest')
clim_mslp_land = clim_mslp.sel(lat=20,lon=80,method='nearest')

clim_skt_ocean = clim_skt.sel(lat=-20,lon=80,method='nearest')
clim_mslp_ocean = clim_mslp.sel(lat=-20,lon=80,method='nearest')


# plotting double y-axis plot with mslp on one side and skin temp on the other

fig, ax1 = plt.subplots(figsize=(16,8))
ax2 = ax1.twinx()
ax1.plot(clim_skt_land-273.5,'b.--',label='Skin Temp')
ax2.plot(clim_mslp_land/100,'r.--',label='MSLP')
ax1.set_xlabel('Time')
ax1.set_ylabel('Skin Temp (degC)')
ax2.set_ylabel('MSLP (hPa)')
ax1.legend(loc='upper left')
ax2.legend(loc='upper right')
plt.title('Time Series of Skin Temp and MSLP at Land (20N,80E)')
plt.savefig('../monsoon/figures/Time Series of Skin Temp and MSLP at Land (20N,80E).png')


# adding a trendline to both the curves and plotting fitting function alongwith the difference of mslp and skin temperature

# fitting a linear function to the data


import numpy as np
from scipy.interpolate import make_interp_spline

# Plotting the time series of skin temp and MSLP for a particular location

clim_skt = xr.open_dataset('../monsoon/datasets/skt.ymonmean.nc')
clim_mslp = xr.open_dataset('../monsoon/datasets/mslp.ymonmean.nc')

clim_mslp = clim_mslp.mslp
clim_skt = clim_skt.skt

# Plotting the time series of skin temp and MSLP for a particular location

clim_skt_land = clim_skt.sel(lat=20,lon=80,method='nearest')
clim_mslp_land = clim_mslp.sel(lat=20,lon=80,method='nearest')

clim_skt_ocean = clim_skt.sel(lat=-20,lon=80,method='nearest')
clim_mslp_ocean = clim_mslp.sel(lat=-20,lon=80,method='nearest')

# create a smoother set of x values
xnew = np.linspace(0, len(clim_skt_land), 300)

# interpolate the skin temperature and MSLP data using a spline function
skt_spline = make_interp_spline(range(len(clim_skt_land)), clim_skt_land-273.5)
mslp_spline = make_interp_spline(range(len(clim_mslp_land)), clim_mslp_land/100)

# calculate the interpolated values using the new set of x values
skt_smooth = skt_spline(xnew)
mslp_smooth = mslp_spline(xnew)

# plot the smoothed curves
fig, ax1 = plt.subplots(figsize=(16,8))
ax2 = ax1.twinx()
ax1.plot(xnew, skt_smooth, 'b-', label='Skin Temp')
ax2.plot(xnew, mslp_smooth, 'r-', label='MSLP')
ax1.set_xlabel('Time')
ax1.set_ylabel('Skin Temp (degC)')
ax2.set_ylabel('MSLP (hPa)')
ax1.legend(loc='upper left')
ax2.legend(loc='upper right')
plt.title('Time Series of Skin Temp and MSLP at Land (20N,80E)')
plt.savefig('../monsoon/figures/Time Series of Skin Temp and MSLP at Land (20N,80E).png')


# Importing required libraries
import numpy as np
import pandas as pd
import xarray as xr
import matplotlib.pyplot as plt
from scipy.interpolate import make_interp_spline, BSpline

# Plotting the time series of skin temp and MSLP for a particular location

clim_skt = xr.open_dataset('../monsoon/datasets/skt.ymonmean.nc')
clim_mslp = xr.open_dataset('../monsoon/datasets/mslp.ymonmean.nc')

clim_mslp = clim_mslp.mslp
clim_skt = clim_skt.skt

# Plotting the time series of skin temp and MSLP for a particular location

clim_skt_land = clim_skt.sel(lat=20,lon=80,method='nearest')
clim_mslp_land = clim_mslp.sel(lat=20,lon=80,method='nearest')

clim_skt_ocean = clim_skt.sel(lat=-20,lon=80,method='nearest')
clim_mslp_ocean = clim_mslp.sel(lat=-20,lon=80,method='nearest')


# plotting double y-axis plot with mslp on one side and skin temp on the other

fig, ax1 = plt.subplots(figsize=(16,8))
ax2 = ax1.twinx()
ax1.plot(clim_skt_land-273.5,'b.--',label='Skin Temp')
ax2.plot(clim_mslp_land/100,'r.--',label='MSLP')

# # plotting the difference of mslp and skin temp
# plt.plot(clim_skt_land+273.5-clim_mslp_land/100 ,'g.--',label='MSLP - Skin Temp')

ax1.set_xlabel('Time')
ax1.set_ylabel('Skin Temp (degC)')
ax2.set_ylabel('MSLP (hPa)')
ax1.legend(loc='upper left')
ax2.legend(loc='upper right')

# Adding smooth curves to the plot

# Define x and y values for skin temp
x1 = np.arange(len(clim_skt_land.time))
y1 = clim_skt_land.values - 273.5

# Define x and y values for MSLP
x2 = np.arange(len(clim_mslp_land.time))
y2 = clim_mslp_land.values / 100

# Smooth curves for skin temp and MSLP
s1 = make_interp_spline(x1, y1, k=3)
s2 = make_interp_spline(x2, y2, k=3)

# Plot the smooth curves
x1_new = np.linspace(x1.min(), x1.max(), 300)
x2_new = np.linspace(x2.min(), x2.max(), 300)

ax1.plot(x1_new, s1(x1_new), 'b-', label='Skin Temp smooth curve')
ax2.plot(x2_new, s2(x2_new), 'r-', label='MSLP smooth curve')

# Add actual data points to the plot
ax1.scatter(x1, y1, s=10, alpha=0.5, color='b', label='Skin Temp data points')
ax2.scatter(x2, y2, s=10, alpha=0.5, color='r', label='MSLP data points')

ax1.legend(loc='center left', frameon=False)
ax2.legend(loc='center right', frameon=False)
plt.title('Time Series of Skin Temp and MSLP at Land (20N,80E)')
plt.savefig('../monsoon/figures/Time Series of Skin Temp and MSLP at Land (20N,80E).png')
plt.show()


# Importing required libraries
import numpy as np
import pandas as pd
import xarray as xr
import matplotlib.pyplot as plt
from scipy.interpolate import make_interp_spline, BSpline

# Plotting the time series of skin temp and MSLP for a particular location

clim_skt = xr.open_dataset('../monsoon/datasets/skt.ymonmean.nc')
clim_mslp = xr.open_dataset('../monsoon/datasets/mslp.ymonmean.nc')

clim_mslp = clim_mslp.mslp
clim_skt = clim_skt.skt

# Plotting the time series of skin temp and MSLP for a particular location

clim_skt_land = clim_skt.sel(lat=20,lon=80,method='nearest')
clim_mslp_land = clim_mslp.sel(lat=20,lon=80,method='nearest')

clim_skt_ocean = clim_skt.sel(lat=-20,lon=80,method='nearest')
clim_mslp_ocean = clim_mslp.sel(lat=-20,lon=80,method='nearest')


# plotting double y-axis plot with mslp on one side and skin temp on the other

fig, ax1 = plt.subplots(figsize=(16,8))
ax2 = ax1.twinx()
ax1.plot(clim_skt_ocean-273.5,'b.--',label='Surface Temperature')
ax2.plot(clim_mslp_ocean/100,'r.--',label='MSLP')


ax1.set_xlabel('Time')
ax1.set_ylabel('Skin Temp (degC)')
ax2.set_ylabel('MSLP (hPa)')
ax1.legend(loc='upper left')
ax2.legend(loc='upper right')

# Adding smooth curves to the plot

# Define x and y values for skin temp
x1 = np.arange(len(clim_skt_ocean.time))
y1 = clim_skt_ocean.values - 273.5

# Define x and y values for MSLP
x2 = np.arange(len(clim_mslp_ocean.time))
y2 = clim_mslp_ocean.values / 100

# Smooth curves for skin temp and MSLP
s1 = make_interp_spline(x1, y1, k=3)
s2 = make_interp_spline(x2, y2, k=3)

# Plot the smooth curves
x1_new = np.linspace(x1.min(), x1.max(), 300)
x2_new = np.linspace(x2.min(), x2.max(), 300)

ax1.plot(x1_new, s1(x1_new), 'b-', label='Skin Temp smooth curve')
ax2.plot(x2_new, s2(x2_new), 'r-', label='MSLP smooth curve')

# Add actual data points to the plot
ax1.scatter(x1, y1, s=10, alpha=0.5, color='b', label='Skin Temp data points')
ax2.scatter(x2, y2, s=10, alpha=0.5, color='r', label='MSLP data points')

ax1.legend(loc='center left', frameon=True)
ax2.legend(loc='upper right', frameon=True)
plt.title('Time Series of Skin Temp and MSLP over Ocean (20S,80E)')
plt.savefig('../monsoon/figures/Time Series of Skin Temp and MSLP over Ocean (20S,80E).png')
plt.show()


month = ['may', 'june', 'july']
level = ['1000','900','850']

# looping through each month and level and plotting the data for each month and level

for i in month:
    for j in level:
        clim_month = xr.open_dataset('../monsoon/datasets/air.'+i+'.nc')
        clim_month = clim_month.mean(dim='time')
        clim_month = clim_month.sel(level=j,method='nearest')
        clim_data, lon = add_cyclic_point(clim_month.air, coord=clim_month.lon)
        lat = clim_month.lat

        fig = plt.figure(figsize=(16,16))
        ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=0.0, globe=None))
        mp = ax.imshow(clim_data-273,extent=(lon.min(),lon.max(),lat.min(),lat.max()),cmap='jet',origin='upper')
        plt.title('Mean air temperature ('+i.capitalize()+')' + ' ' + j + 'hPa')
        ax.add_feature(cfeature.LAND)
        ax.add_feature(cfeature.COASTLINE)
        ax.add_feature(cfeature.OCEAN)
        cbar = fig.colorbar(mp, shrink=0.3,label='Temperature (degC)')
        cbar.minorticks_on()
        gl = ax.gridlines(draw_labels=True,alpha=0.1)
        gl.top_labels = False
        gl.right_labels = False
        plt.savefig('../monsoon/figures/Mean air temperature ('+i+')'+ ' ' +j +'hPa.png')
        plt.show()


from cartopy.util import add_cyclic_point

may = xr.open_dataset('../monsoon/datasets/skt.may.nc')

# may = may.skt
clim_may = may.mean(dim='time')

## adding the cyclic points to remove the missing data at the center of the plot

data = clim_may['skt']
lon = clim_may.coords['lon']

print("Original shape -", data.shape)
lon_idx = data.dims.index('lon')
wrap_data, wrap_lon = add_cyclic_point(data.values, coord=lon, axis=lon_idx)
print("New shape -", wrap_data.shape)


# Plotting the data 

# lon = clim_may.lon
lat = clim_may.lat

fig = plt.figure(figsize=(16,16))
ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=0.0, globe=None))

mp = ax.imshow(wrap_data-273.5,extent=(wrap_lon.min(),wrap_lon.max(),lat.min(),lat.max()),cmap='jet',origin='upper')
plt.title('Long Term Mean Temperature (May)')
# plt.legend(['Temp'])

states_provinces = cfeature.NaturalEarthFeature(
        category='cultural',
        name='admin_1_states_provinces_lines',
        scale='10m',
        facecolor='none')
# ax.add_feature(cfeature.BORDERS,edgecolor='blue')
# ax.add_feature(states_provinces, edgecolor='blue')

ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.OCEAN)


cbar = fig.colorbar(mp, shrink=0.3,label='Temperature (degC)')
cbar.minorticks_on()

#adding the long lat grids and enabling the tick labels
gl = ax.gridlines(draw_labels=True,alpha=0.1)
gl.top_labels = False
gl.right_labels = False

Original shape - (94, 192)
New shape - (94, 193)


uwnd_may =  xr.open_dataset('/home/shiv/Documents/GitHub/EES405/datasets/uwnd.may.850.nc')


uwnd_may = uwnd_may.uwnd


uwnd_may

<xarray.DataArray 'uwnd' (time: 1, level: 1, lat: 73, lon: 144)>
[10512 values with dtype=float32]
Coordinates:
  * time     (time) datetime64[ns] 2000-05-01
  * lon      (lon) float32 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
  * lat      (lat) float32 90.0 87.5 85.0 82.5 80.0 ... -82.5 -85.0 -87.5 -90.0
  * level    (level) float32 850.0
Attributes: (12/13)
    standard_name:  eastward_wind
    long_name:      Monthly U-wind on Pressure Levels
    units:          m/s
    cell_methods:   time: mean (monthly from 6-hourly values)
    precision:      2
    GRIB_id:        33
    ...             ...
    var_desc:       u-wind
    dataset:        NCEP/DOE AMIP-II Reanalysis (Reanalysis-2) Monthly Averages
    level_desc:     Pressure Levels
    statistic:      Individual Obs
    parent_stat:    Other
    actual_range:   [-58.910004 103.159996]

[10512 values with dtype=float32]

array(['2000-05-01T00:00:00.000000000'], dtype='datetime64[ns]')

array([  0. ,   2.5,   5. ,   7.5,  10. ,  12.5,  15. ,  17.5,  20. ,  22.5,
        25. ,  27.5,  30. ,  32.5,  35. ,  37.5,  40. ,  42.5,  45. ,  47.5,
        50. ,  52.5,  55. ,  57.5,  60. ,  62.5,  65. ,  67.5,  70. ,  72.5,
        75. ,  77.5,  80. ,  82.5,  85. ,  87.5,  90. ,  92.5,  95. ,  97.5,
       100. , 102.5, 105. , 107.5, 110. , 112.5, 115. , 117.5, 120. , 122.5,
       125. , 127.5, 130. , 132.5, 135. , 137.5, 140. , 142.5, 145. , 147.5,
       150. , 152.5, 155. , 157.5, 160. , 162.5, 165. , 167.5, 170. , 172.5,
       175. , 177.5, 180. , 182.5, 185. , 187.5, 190. , 192.5, 195. , 197.5,
       200. , 202.5, 205. , 207.5, 210. , 212.5, 215. , 217.5, 220. , 222.5,
       225. , 227.5, 230. , 232.5, 235. , 237.5, 240. , 242.5, 245. , 247.5,
       250. , 252.5, 255. , 257.5, 260. , 262.5, 265. , 267.5, 270. , 272.5,
       275. , 277.5, 280. , 282.5, 285. , 287.5, 290. , 292.5, 295. , 297.5,
       300. , 302.5, 305. , 307.5, 310. , 312.5, 315. , 317.5, 320. , 322.5,
       325. , 327.5, 330. , 332.5, 335. , 337.5, 340. , 342.5, 345. , 347.5,
       350. , 352.5, 355. , 357.5], dtype=float32)

array([ 90. ,  87.5,  85. ,  82.5,  80. ,  77.5,  75. ,  72.5,  70. ,  67.5,
        65. ,  62.5,  60. ,  57.5,  55. ,  52.5,  50. ,  47.5,  45. ,  42.5,
        40. ,  37.5,  35. ,  32.5,  30. ,  27.5,  25. ,  22.5,  20. ,  17.5,
        15. ,  12.5,  10. ,   7.5,   5. ,   2.5,   0. ,  -2.5,  -5. ,  -7.5,
       -10. , -12.5, -15. , -17.5, -20. , -22.5, -25. , -27.5, -30. , -32.5,
       -35. , -37.5, -40. , -42.5, -45. , -47.5, -50. , -52.5, -55. , -57.5,
       -60. , -62.5, -65. , -67.5, -70. , -72.5, -75. , -77.5, -80. , -82.5,
       -85. , -87.5, -90. ], dtype=float32)

array([850.], dtype=float32)

PandasIndex(DatetimeIndex(['2000-05-01'], dtype='datetime64[ns]', name='time', freq=None))

PandasIndex(Float64Index([  0.0,   2.5,   5.0,   7.5,  10.0,  12.5,  15.0,  17.5,  20.0,
               22.5,
              ...
              335.0, 337.5, 340.0, 342.5, 345.0, 347.5, 350.0, 352.5, 355.0,
              357.5],
             dtype='float64', name='lon', length=144))

PandasIndex(Float64Index([ 90.0,  87.5,  85.0,  82.5,  80.0,  77.5,  75.0,  72.5,  70.0,
               67.5,  65.0,  62.5,  60.0,  57.5,  55.0,  52.5,  50.0,  47.5,
               45.0,  42.5,  40.0,  37.5,  35.0,  32.5,  30.0,  27.5,  25.0,
               22.5,  20.0,  17.5,  15.0,  12.5,  10.0,   7.5,   5.0,   2.5,
                0.0,  -2.5,  -5.0,  -7.5, -10.0, -12.5, -15.0, -17.5, -20.0,
              -22.5, -25.0, -27.5, -30.0, -32.5, -35.0, -37.5, -40.0, -42.5,
              -45.0, -47.5, -50.0, -52.5, -55.0, -57.5, -60.0, -62.5, -65.0,
              -67.5, -70.0, -72.5, -75.0, -77.5, -80.0, -82.5, -85.0, -87.5,
              -90.0],
             dtype='float64', name='lat'))

PandasIndex(Float64Index([850.0], dtype='float64', name='level'))


uwnd_may = uwnd_may.squeeze(dim=['time','level'])


uwnd_may

<xarray.DataArray 'uwnd' (lat: 73, lon: 144)>
[10512 values with dtype=float32]
Coordinates:
    time     datetime64[ns] 2000-05-01
  * lon      (lon) float32 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
  * lat      (lat) float32 90.0 87.5 85.0 82.5 80.0 ... -82.5 -85.0 -87.5 -90.0
    level    float32 850.0
Attributes: (12/13)
    standard_name:  eastward_wind
    long_name:      Monthly U-wind on Pressure Levels
    units:          m/s
    cell_methods:   time: mean (monthly from 6-hourly values)
    precision:      2
    GRIB_id:        33
    ...             ...
    var_desc:       u-wind
    dataset:        NCEP/DOE AMIP-II Reanalysis (Reanalysis-2) Monthly Averages
    level_desc:     Pressure Levels
    statistic:      Individual Obs
    parent_stat:    Other
    actual_range:   [-58.910004 103.159996]

[10512 values with dtype=float32]

array('2000-05-01T00:00:00.000000000', dtype='datetime64[ns]')

array([  0. ,   2.5,   5. ,   7.5,  10. ,  12.5,  15. ,  17.5,  20. ,  22.5,
        25. ,  27.5,  30. ,  32.5,  35. ,  37.5,  40. ,  42.5,  45. ,  47.5,
        50. ,  52.5,  55. ,  57.5,  60. ,  62.5,  65. ,  67.5,  70. ,  72.5,
        75. ,  77.5,  80. ,  82.5,  85. ,  87.5,  90. ,  92.5,  95. ,  97.5,
       100. , 102.5, 105. , 107.5, 110. , 112.5, 115. , 117.5, 120. , 122.5,
       125. , 127.5, 130. , 132.5, 135. , 137.5, 140. , 142.5, 145. , 147.5,
       150. , 152.5, 155. , 157.5, 160. , 162.5, 165. , 167.5, 170. , 172.5,
       175. , 177.5, 180. , 182.5, 185. , 187.5, 190. , 192.5, 195. , 197.5,
       200. , 202.5, 205. , 207.5, 210. , 212.5, 215. , 217.5, 220. , 222.5,
       225. , 227.5, 230. , 232.5, 235. , 237.5, 240. , 242.5, 245. , 247.5,
       250. , 252.5, 255. , 257.5, 260. , 262.5, 265. , 267.5, 270. , 272.5,
       275. , 277.5, 280. , 282.5, 285. , 287.5, 290. , 292.5, 295. , 297.5,
       300. , 302.5, 305. , 307.5, 310. , 312.5, 315. , 317.5, 320. , 322.5,
       325. , 327.5, 330. , 332.5, 335. , 337.5, 340. , 342.5, 345. , 347.5,
       350. , 352.5, 355. , 357.5], dtype=float32)

array([ 90. ,  87.5,  85. ,  82.5,  80. ,  77.5,  75. ,  72.5,  70. ,  67.5,
        65. ,  62.5,  60. ,  57.5,  55. ,  52.5,  50. ,  47.5,  45. ,  42.5,
        40. ,  37.5,  35. ,  32.5,  30. ,  27.5,  25. ,  22.5,  20. ,  17.5,
        15. ,  12.5,  10. ,   7.5,   5. ,   2.5,   0. ,  -2.5,  -5. ,  -7.5,
       -10. , -12.5, -15. , -17.5, -20. , -22.5, -25. , -27.5, -30. , -32.5,
       -35. , -37.5, -40. , -42.5, -45. , -47.5, -50. , -52.5, -55. , -57.5,
       -60. , -62.5, -65. , -67.5, -70. , -72.5, -75. , -77.5, -80. , -82.5,
       -85. , -87.5, -90. ], dtype=float32)

array(850., dtype=float32)

PandasIndex(Float64Index([  0.0,   2.5,   5.0,   7.5,  10.0,  12.5,  15.0,  17.5,  20.0,
               22.5,
              ...
              335.0, 337.5, 340.0, 342.5, 345.0, 347.5, 350.0, 352.5, 355.0,
              357.5],
             dtype='float64', name='lon', length=144))

PandasIndex(Float64Index([ 90.0,  87.5,  85.0,  82.5,  80.0,  77.5,  75.0,  72.5,  70.0,
               67.5,  65.0,  62.5,  60.0,  57.5,  55.0,  52.5,  50.0,  47.5,
               45.0,  42.5,  40.0,  37.5,  35.0,  32.5,  30.0,  27.5,  25.0,
               22.5,  20.0,  17.5,  15.0,  12.5,  10.0,   7.5,   5.0,   2.5,
                0.0,  -2.5,  -5.0,  -7.5, -10.0, -12.5, -15.0, -17.5, -20.0,
              -22.5, -25.0, -27.5, -30.0, -32.5, -35.0, -37.5, -40.0, -42.5,
              -45.0, -47.5, -50.0, -52.5, -55.0, -57.5, -60.0, -62.5, -65.0,
              -67.5, -70.0, -72.5, -75.0, -77.5, -80.0, -82.5, -85.0, -87.5,
              -90.0],
             dtype='float64', name='lat'))


uwnd_may.shape

(73, 144)


from cartopy.util import add_cyclic_point

# may = xr.open_dataset('../monsoon/datasets/skt.may.nc')

# may = may.skt
# clim_may = may.mean(dim='time')

## adding the cyclic points to remove the missing data at the center of the plot

# data = clim_may['uwnd']
lon = uwnd_may.coords['lon']

print("Original shape -", uwnd_may.shape)
lon_idx = data.dims.index('lon')
wrap_data_uwnd_may, wrap_lon = add_cyclic_point(uwnd_may.values, coord=lon, axis=lon_idx)
print("New shape -", wrap_data.shape)


# Plotting the data 

# lon = uwnd_may.lon
lat = uwnd_may.lat

fig = plt.figure(figsize=(16,16))
ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=0.0, globe=None))

mp = ax.imshow(wrap_data_uwnd_may,extent=(wrap_lon.min(),wrap_lon.max(),lat.min(),lat.max()),cmap='jet',origin='upper')
plt.title('Long Term Mean Temperature (May)')
# plt.legend(['Temp'])

states_provinces = cfeature.NaturalEarthFeature(
        category='cultural',
        name='admin_1_states_provinces_lines',
        scale='10m',
        facecolor='none')
# ax.add_feature(cfeature.BORDERS,edgecolor='blue')
# ax.add_feature(states_provinces, edgecolor='blue')

ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.OCEAN)


cbar = fig.colorbar(mp, shrink=0.3,label='Wind Speed')
cbar.minorticks_on()

#adding the long lat grids and enabling the tick labels
gl = ax.gridlines(draw_labels=True,alpha=0.1)
gl.top_labels = False
gl.right_labels = False

Original shape - (73, 144)
New shape - (94, 193)


vwnd_may = xr.open_dataset('/home/shiv/Documents/GitHub/EES405/datasets/vwnd.may.850.nc')
vwnd_may = vwnd_may.vwnd
vwnd_may = vwnd_may.squeeze(dim=['time','level'])


from cartopy.util import add_cyclic_point

# may = xr.open_dataset('../monsoon/datasets/skt.may.nc')

# may = may.skt
# clim_may = may.mean(dim='time')

## adding the cyclic points to remove the missing data at the center of the plot

# data = clim_may['uwnd']
lon = uwnd_may.coords['lon']

print("Original shape -", uwnd_may.shape)

# Add cyclic point to uwnd and vwnd
lon_idx = uwnd_may.dims.index('lon')
wrap_data_uwnd_may, wrap_lon = add_cyclic_point(uwnd_may.values, coord=lon, axis=lon_idx)
wrap_data_vwnd_may,wrap_lon = add_cyclic_point(vwnd_may.values, coord=lon, axis=lon_idx)

# Compute vector addition
wind_may = np.sqrt(wrap_data_uwnd_may**2 + wrap_data_vwnd_may**2)
wind_may_direction = np.arctan2(wrap_data_vwnd_may, wrap_data_uwnd_may)

u_wind_add = wind_may * np.cos(wind_may_direction)
v_wind_add = wind_may * np.sin(wind_may_direction)

print("New shape -", wrap_data_uwnd_may.shape)

# Plotting the data 

# lon = uwnd_may.lon
lat = uwnd_may.lat

fig = plt.figure(figsize=(16,16))
ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=0.0, globe=None))

# Create 2D grid of longitudes and latitudes
lon_2d, lat_2d = np.meshgrid(wrap_lon, lat)

# Compute vector addition of u-wind and v-wind
u_wind_add = u_wind_add + v_wind_add
v_wind_add = v_wind_add + u_wind_add

# Compute wind speed from u-wind and v-wind
wind_may = np.sqrt(u_wind_add**2 + v_wind_add**2)


# Plot vector addition
skip = (slice(None, None, 1), slice(None, None, 4))
ax.quiver(lon_2d[skip], lat_2d[skip], u_wind_add[skip], v_wind_add[skip], wind_may[skip],
          transform=ccrs.PlateCarree(), cmap='coolwarm', scale=300)



# Plot wind speed as contour plot
mp = ax.imshow(wind_may, extent=(wrap_lon.min(),wrap_lon.max(),lat.min(),lat.max()), cmap='jet', origin='upper')


states_provinces = cfeature.NaturalEarthFeature(
        category='cultural',
        name='admin_1_states_provinces_lines',
        scale='10m',
        facecolor='none')
# ax.add_feature(cfeature.BORDERS,edgecolor='blue')
# ax.add_feature(states_provinces, edgecolor='blue')

ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.OCEAN)


cbar = fig.colorbar(mp, shrink=0.3,label='Wind Speed')
cbar.minorticks_on()

#adding the long lat grids and enabling the tick labels
gl = ax.gridlines(draw_labels=True,alpha=0.1)
gl.top_labels = False
gl.right_labels = False

Original shape - (73, 144)
New shape - (73, 145)


# making all the plots at once for may, june, july, august, september and jjas 

months = ['may','june','july','august','september','jjas']

for i in months:
    uwnd = xr.open_dataset('/home/shiv/Documents/GitHub/EES405/datasets/uwnd.'+i+'.850.nc')
    vwnd = xr.open_dataset('/home/shiv/Documents/GitHub/EES405/datasets/vwnd.'+i+'.850.nc')
    uwnd = uwnd.uwnd
    vwnd = vwnd.vwnd
    uwnd = uwnd.squeeze(dim=['time','level'])
    vwnd = vwnd.squeeze(dim=['time','level'])
    
    from cartopy.util import add_cyclic_point

    lon = uwnd_may.coords['lon']

    print("Original shape -", uwnd_may.shape)

    # Add cyclic point to uwnd and vwnd
    lon_idx = uwnd_may.dims.index('lon')
    wrap_data_uwnd_may, wrap_lon = add_cyclic_point(uwnd_may.values, coord=lon, axis=lon_idx)
    wrap_data_vwnd_may,wrap_lon = add_cyclic_point(vwnd_may.values, coord=lon, axis=lon_idx)

    # Compute vector addition
    wind_may = np.sqrt(wrap_data_uwnd_may**2 + wrap_data_vwnd_may**2)
    wind_may_direction = np.arctan2(wrap_data_vwnd_may, wrap_data_uwnd_may)

    u_wind_add = wind_may * np.cos(wind_may_direction)
    v_wind_add = wind_may * np.sin(wind_may_direction)

    print("New shape -", wrap_data_uwnd_may.shape)

    # Plotting the data 

    # lon = uwnd_may.lon
    lat = uwnd_may.lat

    fig = plt.figure(figsize=(16,16))
    ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=0.0, globe=None))

    # Create 2D grid of longitudes and latitudes
    lon_2d, lat_2d = np.meshgrid(wrap_lon, lat)

    # Compute vector addition of u-wind and v-wind
    u_wind_add = u_wind_add + v_wind_add
    v_wind_add = v_wind_add + u_wind_add

    # Compute wind speed from u-wind and v-wind
    wind_may = np.sqrt(u_wind_add**2 + v_wind_add**2)


    # Plot vector addition
    skip = (slice(None, None, 1), slice(None, None, 4))
    ax.quiver(lon_2d[skip], lat_2d[skip], u_wind_add[skip], v_wind_add[skip], wind_may[skip],
            transform=ccrs.PlateCarree(), cmap='coolwarm', scale=300)



    # Plot wind speed as contour plot
    mp = ax.imshow(wind_may, extent=(wrap_lon.min(),wrap_lon.max(),lat.min(),lat.max()), cmap='jet', origin='upper')
    plt.title('Wind Speed and Direction for '+i + ' at 850 hPa')



    states_provinces = cfeature.NaturalEarthFeature(
            category='cultural',
            name='admin_1_states_provinces_lines',
            scale='10m',
            facecolor='none')
    # ax.add_feature(cfeature.BORDERS,edgecolor='blue')
    # ax.add_feature(states_provinces, edgecolor='blue')

    ax.add_feature(cfeature.LAND)
    ax.add_feature(cfeature.COASTLINE)
    ax.add_feature(cfeature.OCEAN)


    cbar = fig.colorbar(mp, shrink=0.3,label='Wind Speed')
    cbar.minorticks_on()

    #adding the long lat grids and enabling the tick labels
    gl = ax.gridlines(draw_labels=True,alpha=0.1)
    gl.top_labels = False
    gl.right_labels = False

    plt.savefig('/home/shiv/Documents/GitHub/EES405/monsoon/figures/'+i+'850hPa.png',dpi=300)

Original shape - (73, 144)
New shape - (73, 145)
Original shape - (73, 144)
New shape - (73, 145)
Original shape - (73, 144)
New shape - (73, 145)
Original shape - (73, 144)
New shape - (73, 145)
Original shape - (73, 144)
New shape - (73, 145)
Original shape - (73, 144)
New shape - (73, 145)


uwnd_data = xr.open_dataset('/home/shiv/Documents/GitHub/EES405/datasets/uwnd.may.nc')
vwnd_data = xr.open_dataset('/home/shiv/Documents/GitHub/EES405/datasets/vwnd.may.nc')

uwnd_data = uwnd_data.uwnd
vwnd_data = vwnd_data.vwnd
uwnd_data = uwnd_data.sel(level=200)
vwnd_data = vwnd_data.sel(level=200)
uwnd_data = uwnd_data.squeeze(dim=['time'])
vwnd_data = vwnd_data.squeeze(dim=['time'])


vwnd_data = vwnd_data.sel(lat=slice(0,90))

<xarray.DataArray 'vwnd' (lat: 73, lon: 144)>
[10512 values with dtype=float32]
Coordinates:
    time     datetime64[ns] 2000-05-01
  * lon      (lon) float32 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
  * lat      (lat) float32 90.0 87.5 85.0 82.5 80.0 ... -82.5 -85.0 -87.5 -90.0
    level    float32 200.0
Attributes: (12/13)
    standard_name:  northward_wind
    long_name:      Monthly V-wind on Pressure Levels
    units:          m/s
    cell_methods:   time: mean (monthly from 6-hourly values) Monthly Averages
    precision:      2
    GRIB_id:        34
    ...             ...
    var_desc:       v-wind
    dataset:        NCEP/DOE AMIP-II Reanalysis (Reanalysis-2) Monthly Averages
    level_desc:     Pressure Levels
    statistic:      Individual Obs
    parent_stat:    Other
    actual_range:   [-69.100006  69.21    ]

[10512 values with dtype=float32]

array('2000-05-01T00:00:00.000000000', dtype='datetime64[ns]')

array([  0. ,   2.5,   5. ,   7.5,  10. ,  12.5,  15. ,  17.5,  20. ,  22.5,
        25. ,  27.5,  30. ,  32.5,  35. ,  37.5,  40. ,  42.5,  45. ,  47.5,
        50. ,  52.5,  55. ,  57.5,  60. ,  62.5,  65. ,  67.5,  70. ,  72.5,
        75. ,  77.5,  80. ,  82.5,  85. ,  87.5,  90. ,  92.5,  95. ,  97.5,
       100. , 102.5, 105. , 107.5, 110. , 112.5, 115. , 117.5, 120. , 122.5,
       125. , 127.5, 130. , 132.5, 135. , 137.5, 140. , 142.5, 145. , 147.5,
       150. , 152.5, 155. , 157.5, 160. , 162.5, 165. , 167.5, 170. , 172.5,
       175. , 177.5, 180. , 182.5, 185. , 187.5, 190. , 192.5, 195. , 197.5,
       200. , 202.5, 205. , 207.5, 210. , 212.5, 215. , 217.5, 220. , 222.5,
       225. , 227.5, 230. , 232.5, 235. , 237.5, 240. , 242.5, 245. , 247.5,
       250. , 252.5, 255. , 257.5, 260. , 262.5, 265. , 267.5, 270. , 272.5,
       275. , 277.5, 280. , 282.5, 285. , 287.5, 290. , 292.5, 295. , 297.5,
       300. , 302.5, 305. , 307.5, 310. , 312.5, 315. , 317.5, 320. , 322.5,
       325. , 327.5, 330. , 332.5, 335. , 337.5, 340. , 342.5, 345. , 347.5,
       350. , 352.5, 355. , 357.5], dtype=float32)

array([ 90. ,  87.5,  85. ,  82.5,  80. ,  77.5,  75. ,  72.5,  70. ,  67.5,
        65. ,  62.5,  60. ,  57.5,  55. ,  52.5,  50. ,  47.5,  45. ,  42.5,
        40. ,  37.5,  35. ,  32.5,  30. ,  27.5,  25. ,  22.5,  20. ,  17.5,
        15. ,  12.5,  10. ,   7.5,   5. ,   2.5,   0. ,  -2.5,  -5. ,  -7.5,
       -10. , -12.5, -15. , -17.5, -20. , -22.5, -25. , -27.5, -30. , -32.5,
       -35. , -37.5, -40. , -42.5, -45. , -47.5, -50. , -52.5, -55. , -57.5,
       -60. , -62.5, -65. , -67.5, -70. , -72.5, -75. , -77.5, -80. , -82.5,
       -85. , -87.5, -90. ], dtype=float32)

array(200., dtype=float32)

PandasIndex(Float64Index([  0.0,   2.5,   5.0,   7.5,  10.0,  12.5,  15.0,  17.5,  20.0,
               22.5,
              ...
              335.0, 337.5, 340.0, 342.5, 345.0, 347.5, 350.0, 352.5, 355.0,
              357.5],
             dtype='float64', name='lon', length=144))

PandasIndex(Float64Index([ 90.0,  87.5,  85.0,  82.5,  80.0,  77.5,  75.0,  72.5,  70.0,
               67.5,  65.0,  62.5,  60.0,  57.5,  55.0,  52.5,  50.0,  47.5,
               45.0,  42.5,  40.0,  37.5,  35.0,  32.5,  30.0,  27.5,  25.0,
               22.5,  20.0,  17.5,  15.0,  12.5,  10.0,   7.5,   5.0,   2.5,
                0.0,  -2.5,  -5.0,  -7.5, -10.0, -12.5, -15.0, -17.5, -20.0,
              -22.5, -25.0, -27.5, -30.0, -32.5, -35.0, -37.5, -40.0, -42.5,
              -45.0, -47.5, -50.0, -52.5, -55.0, -57.5, -60.0, -62.5, -65.0,
              -67.5, -70.0, -72.5, -75.0, -77.5, -80.0, -82.5, -85.0, -87.5,
              -90.0],
             dtype='float64', name='lat'))


uwnd_data = xr.open_dataset('/home/shiv/Documents/GitHub/EES405/datasets/uwnd.may.nc')
vwnd_data = xr.open_dataset('/home/shiv/Documents/GitHub/EES405/datasets/vwnd.may.nc')

uwnd_data = uwnd_data.uwnd
vwnd_data = vwnd_data.vwnd
uwnd_data = uwnd_data.sel(level=200)
vwnd_data = vwnd_data.sel(level=200)
uwnd_data = uwnd_data.squeeze(dim=['time'])
vwnd_data = vwnd_data.squeeze(dim=['time'])

lon = uwnd_data.coords['lon']

print("Original shape -", uwnd_data.shape)

# Add cyclic point to uwnd and vwnd
lon_idx = uwnd_data.dims.index('lon')
wrap_data_uwnd_data, wrap_lon = add_cyclic_point(uwnd_data.values, coord=lon, axis=lon_idx)
wrap_data_vwnd_data, wrap_lon = add_cyclic_point(vwnd_data.values, coord=lon, axis=lon_idx)

# Compute vector addition
wind_may = np.sqrt(wrap_data_uwnd_data**2 + wrap_data_vwnd_data**2)
wind_may_direction = np.arctan2(wrap_data_vwnd_data, wrap_data_uwnd_data)

u_wind_add = wind_may * np.cos(wind_may_direction)
v_wind_add = wind_may * np.sin(wind_may_direction)

print("New shape -", wrap_data_uwnd_may.shape)

# Plotting the data 

# lon = uwnd_may.lon
lat = uwnd_may.lat

fig = plt.figure(figsize=(16,16))
ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=0.0, globe=None))

# Create 2D grid of longitudes and latitudes
lon_2d, lat_2d = np.meshgrid(wrap_lon, lat)

# Compute vector addition of u-wind and v-wind
u_wind_add = u_wind_add + v_wind_add
v_wind_add = v_wind_add + u_wind_add

# Compute wind speed from u-wind and v-wind
wind_may = np.sqrt(u_wind_add**2 + v_wind_add**2)


# Plot vector addition
skip = (slice(None, None, 3), slice(None, None, 4))
ax.quiver(lon_2d[skip], lat_2d[skip], u_wind_add[skip], v_wind_add[skip], wind_may[skip],
          transform=ccrs.PlateCarree(), cmap='coolwarm', scale=500)



# Plot wind speed as contour plot
mp = ax.imshow(wind_may, extent=(wrap_lon.min(),wrap_lon.max(),lat.min(),lat.max()), cmap='jet', origin='upper')


states_provinces = cfeature.NaturalEarthFeature(
        category='cultural',
        name='admin_1_states_provinces_lines',
        scale='10m',
        facecolor='none')
# ax.add_feature(cfeature.BORDERS,edgecolor='blue')
# ax.add_feature(states_provinces, edgecolor='blue')

ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.OCEAN)


cbar = fig.colorbar(mp, shrink=0.3,label='Wind Speed')
cbar.minorticks_on()

#adding the long lat grids and enabling the tick labels
gl = ax.gridlines(draw_labels=True,alpha=0.1)
gl.top_labels = False
gl.right_labels = False

Original shape - (73, 144)
New shape - (73, 145)


uwnd_data = xr.open_dataset('/home/shiv/Documents/GitHub/EES405/datasets/uwnd.may.nc')
vwnd_data = xr.open_dataset('/home/shiv/Documents/GitHub/EES405/datasets/vwnd.may.nc')

uwnd_data = uwnd_data.uwnd
vwnd_data = vwnd_data.vwnd
uwnd_data = uwnd_data.sel(level=200)
vwnd_data = vwnd_data.sel(level=200)
uwnd_data = uwnd_data.squeeze(dim=['time'])
vwnd_data = vwnd_data.squeeze(dim=['time'])

lon = uwnd_data.coords['lon']

print("Original shape -", uwnd_data.shape)

# Add cyclic point to uwnd and vwnd
lon_idx = uwnd_data.dims.index('lon')
wrap_data_uwnd_data, wrap_lon = add_cyclic_point(uwnd_data.values, coord=lon, axis=lon_idx)
wrap_data_vwnd_data, wrap_lon = add_cyclic_point(vwnd_data.values, coord=lon, axis=lon_idx)

# Compute vector addition
wind_may = np.sqrt(wrap_data_uwnd_data**2 + wrap_data_vwnd_data**2)
wind_may_direction = np.arctan2(wrap_data_vwnd_data, wrap_data_uwnd_data)

u_wind_add = wind_may * np.cos(wind_may_direction)
v_wind_add = wind_may * np.sin(wind_may_direction)

print("New shape -", wrap_data_uwnd_may.shape)

# Plotting the data 

# lon = uwnd_may.lon
lat = uwnd_may.lat

fig = plt.figure(figsize=(16,16))
ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=0.0, globe=None))

# Create 2D grid of longitudes and latitudes
lon_2d, lat_2d = np.meshgrid(wrap_lon, lat)

# Compute vector addition of u-wind and v-wind
u_wind_add = u_wind_add + v_wind_add
v_wind_add = v_wind_add + u_wind_add

# Compute wind speed from u-wind and v-wind
wind_may = np.sqrt(u_wind_add**2 + v_wind_add**2)


# Find the indices of the maximum wind speed
max_idx = np.unravel_index(np.argmax(wind_may), wind_may.shape)

# Plot a red dot at the location of maximum wind speed
ax.plot(lon_2d[max_idx], lat_2d[max_idx], 'ko', markersize=10, transform=ccrs.PlateCarree())

# Plot vector addition
skip = (slice(None, None, 3), slice(None, None, 4))
q = ax.quiver(lon_2d[skip], lat_2d[skip], u_wind_add[skip], v_wind_add[skip], wind_may[skip],
          transform=ccrs.PlateCarree(), cmap='coolwarm', scale=500)

# Plot wind speed as contour plot
mp = ax.imshow(wind_may, extent=(wrap_lon.min(),wrap_lon.max(),lat.min(),lat.max()), cmap='jet', origin='upper')

# Add colorbar
cbar = fig.colorbar(mp, shrink=0.3,label='Wind Speed')
cbar.minorticks_on()

# Add title
ax.set_title('May Wind Speed at 200 hPa')

#adding the long lat grids and enabling the tick labels
gl = ax.gridlines(draw_labels=True,alpha=0.1)
gl.top_labels = False
gl.right_labels = False

# Add land, coastlines, and ocean features
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.OCEAN)
plt.show()

Original shape - (73, 144)
New shape - (73, 145)


# plot all the months together in one code 

months = ['may','june','july','august','september','jjas']

for i in months:
    uwnd_data = xr.open_dataset('/home/shiv/Documents/GitHub/EES405/datasets/uwnd.'+i+'.nc')
    vwnd_data = xr.open_dataset('/home/shiv/Documents/GitHub/EES405/datasets/vwnd.'+i+'.nc')

    uwnd_data = uwnd_data.uwnd
    vwnd_data = vwnd_data.vwnd
    uwnd_data = uwnd_data.sel(level=200)
    vwnd_data = vwnd_data.sel(level=200)
    uwnd_data = uwnd_data.squeeze(dim=['time'])
    vwnd_data = vwnd_data.squeeze(dim=['time'])

    lon = uwnd_data.coords['lon']

    print("Original shape -", uwnd_data.shape)

    # Add cyclic point to uwnd and vwnd
    lon_idx = uwnd_data.dims.index('lon')
    wrap_data_uwnd_data, wrap_lon = add_cyclic_point(uwnd_data.values, coord=lon, axis=lon_idx)
    wrap_data_vwnd_data, wrap_lon = add_cyclic_point(vwnd_data.values, coord=lon, axis=lon_idx)

    # Compute vector addition
    wind_may = np.sqrt(wrap_data_uwnd_data**2 + wrap_data_vwnd_data**2)
    wind_may_direction = np.arctan2(wrap_data_vwnd_data, wrap_data_uwnd_data)

    u_wind_add = wind_may * np.cos(wind_may_direction)
    v_wind_add = wind_may * np.sin(wind_may_direction)

    print("New shape -", wrap_data_uwnd_may.shape)

    # Plotting the data 

    # lon = uwnd_may.lon
    lat = uwnd_may.lat

    fig = plt.figure(figsize=(16,16))
    ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=0.0, globe=None))

    # Create 2D grid of longitudes and latitudes
    lon_2d, lat_2d = np.meshgrid(wrap_lon, lat)

    # Compute vector addition of u-wind and v-wind
    u_wind_add = u_wind_add + v_wind_add
    v_wind_add = v_wind_add + u_wind_add

    # Compute vector addition of u-wind and v-wind
    u_wind_add = u_wind_add + v_wind_add
    v_wind_add = v_wind_add + u_wind_add

    # Compute wind speed from u-wind and v-wind
    wind_may = np.sqrt(u_wind_add**2 + v_wind_add**2)


    # Find the indices of the maximum wind speed
    max_idx = np.unravel_index(np.argmax(wind_may), wind_may.shape)

    # Plot a red dot at the location of maximum wind speed
    ax.plot(lon_2d[max_idx], lat_2d[max_idx], 'ko', markersize=10, transform=ccrs.PlateCarree())

    # Plot vector addition
    skip = (slice(None, None, 3), slice(None, None, 4))
    q = ax.quiver(lon_2d[skip], lat_2d[skip], u_wind_add[skip], v_wind_add[skip], wind_may[skip],
            transform=ccrs.PlateCarree(), cmap='coolwarm', scale=500)

    # Plot wind speed as contour plot
    mp = ax.imshow(wind_may, extent=(wrap_lon.min(),wrap_lon.max(),lat.min(),lat.max()), cmap='jet', origin='upper')

    # Add colorbar
    cbar = fig.colorbar(mp, shrink=0.3,label='Wind Speed')
    cbar.minorticks_on()

    # Add title
    ax.set_title(i +' Wind Speed at 200 hPa')

    #adding the long lat grids and enabling the tick labels
    gl = ax.gridlines(draw_labels=True,alpha=0.1)
    gl.top_labels = False
    gl.right_labels = False

    # Add land, coastlines, and ocean features
    ax.add_feature(cfeature.LAND)
    ax.add_feature(cfeature.COASTLINE)
    ax.add_feature(cfeature.OCEAN)
    plt.show()

Original shape - (73, 144)
New shape - (73, 145)

Original shape - (73, 144)
New shape - (73, 145)

Original shape - (73, 144)
New shape - (73, 145)

Original shape - (73, 144)
New shape - (73, 145)

Original shape - (73, 144)
New shape - (73, 145)

Original shape - (73, 144)
New shape - (73, 145)